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Abstract 

In this paper we consider the suppression of bottomonium states in ultrarelativistic heavy 
ion collisions. We compute the suppression as a function of centrality, rapidity, and trans- 
verse momentum for the states T(ls), T(2s), T(3s), %bii an d X&2- Using this information, 
we then compute the inclusive T(ls) suppression as a function of centrality, rapidity, and 
transverse momentum including feed down effects. Calculations are performed for both 
RHIC = 200 GeV Au-Au collisions and LHC v /s^7 = 2.76 TeV Pb-Pb collisions. 

From the comparison of our theoretical results with data available from the STAR and CMS 
Collaborations we are able to constrain the shear viscosity to entropy ratio to be in the 
range 0.08 < rj/S < 0.24. Our results are consistent with the creation of a high temperature 
quark-gluon plasma at both RHIC and LHC collision energies. 

Keywords: Quarkonium Suppression, Bottomonium Suppression, Relativistic Heavy Ion 
Collision, Quark-Gluon Plasma 



1. Introduction 

The goal of ultrarelativistic heavy ion collision experiments at the Relativistic Heavy 
Ion Collider at Brookhaven National Laboratory (RHIC) and the Large Hadron Collider 
(LHC) at CERN is to create a tiny volume of matter (~ 1000 fm 3 ) which has been heated 
to a temperature exceeding that necessary to deconfine quarks and gluons. Lattice quantum 
chromodynamics (lattice QCD) measurements of the equation of state of strongly interacting 
matter [THE] show that there is crossover from hadronic matter to a quark-gluon plasma at 
temperatures on the order of 175 MeV which corresponds to approximately two trillion 
degrees Kelvin. For RHIC ^s^n — 200 GeV Au-Au collisions, initial maximum central 
temperatures of T ~ 450 MeV were generated and for current LHC ^/s~N r N = 2.76 TeV 
collisions one obtains T ~ 550 MeV [5]. For the upcoming full energy LHC heavy ion runs 
with ^fs^ = 5.5 TeV one expects T ~ 700 - 800 MeV. 

At such extremely high temperatures strongly interacting matter makes a phase transi- 
tion to a deconfmed plasma of quarks and gluons and, as a result, one expects the emergence 
of Debye screening of the interaction between quarks and gluons. This leads to the dissolu- 
tion of hadronic bound states [7j. A particularly interesting subset of hadronic states consists 
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of those which are comprised of heavy quarks because the spectrum of such states can be 
found using potential-based non-relativistic treatments. Based on such potential models 
there were early predictions [HI E] that J/ty production would be suppressed in heavy ion 
collisions relative to the corresponding production in proton-proton collisions scaled by the 
number of nucleons participating in the collision. 

As mentioned above, heavy quarkonium has received the most theoretical attention, 
since heavy quark states are dominated by short distance physics and can be treated using 
heavy quark effective theory. Based on such effective theories of QCD, non-relativistic 
quarkonium states can be reliably described. Their binding energies are much smaller than 
the quark mass tuq 3> Aqcd (Q = c,b), and their sizes are much larger than 1/mg. At 
zero temperature, since the velocity of the quarks in the bound state is small (v <C c), 
quarkonium can be understood in terms of non-relativistic potential models such as the 
Cornell potential which can be derived directly from QCD using effective field theory [T0HT2"] . 
Using such non-relativistic potential models studies of quarkonium spectral functions and 
meson current correlators have been performed [T314T9] . The results have been compared to 
first-principles lattice QCD calculations [201426] which rely on the maximum entropy method 
[2"7l [28J. Additionally, there have been some lattice developments using non-relativistic 
lattice QCD [29]. 

Additionally, in recent years there has been an important theoretical development, 
namely the first-principles calculation of the thermal widths of heavy quarkonium states 
which emerge from imaginary-valued contributions to the heavy quark potential. The first 
calculation of the leading-order perturbative imaginary part of the potential due to glu- 
onic Landau damping was performed by Laine et al. [3"0"1 I3"T] . Since then an additional 
imaginary-valued contribution to the potential coming from singlet to octet transitions has 
also been computed using the effective field theory approach [32], and lattice calculations 
have been performed in order to determine the imaginary part of the heavy quark poten- 
tial [33J. These imaginary contributions to the potential are related to quarkonium decay 
processes in the plasma. The consequences of such imaginary parts on heavy quarkonium 
spectral functions [SUES], perturbative thermal widths [301 136], quarkonium suppression in 
a T-matrix approach [3TH39] . and in stochastic real-time dynamics [40 J have recently been 
studied; however, these studies were restricted to the case of an isotropicthermal plasma, 
which is only the case if one assumes ideal hydrodynamical evolution. 

The calculation of the heavy quark potential has since been extended to the case of a 
plasma with finite momentum-space anisotropy. Both the real [4T1 142] and imaginary [43TI45] 
parts have been computed in this case. Additionally, the impact of the imaginary part of the 
potential on the thermal widths of the states in both isotropic and anisotropic plasmas was 
recently studied [16]. The consideration of momentum-space anisotropic plasmas is necessary 
since, for any finite shear viscosity, the quark-gluon plasma possesses local momentum-space 
anisotropies [47H55] • Depending on the magnitude of the shear viscosity, these momentum- 
space anisotropies can persist for a long time and can be quite large, particularly at early 
times or near the edges of the plasma. This is true for both strong and weak coupling values 
of the shear viscosity and the magnitude of the maximal momentum-space anisotropies 
increases with increasing shear viscosity. In fact, the magnitude of these momentum space 
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anisotropics can become so large that they call into doubt the reliability of the viscous 
hydrodynamical treatment which implicitly assumes a nearly isotropic state. 

This has motivated the development of a new dynamical formalism called "anisotropic 
hydrodynamics" (aHydro) which extends traditional viscous hydrodynamical treatments 
to cases in which the local momentum-space anisotropy of the plasma can be large [511 - 
153] . The result is a dynamical framework that reduces to 2nd order viscous hydrodynamics 
for weakly anisotropic plasmas, but can better describe highly anisotropic plasmas. For 
one-dimensional dynamics which is homogeneous in the transverse directions, the aHydro 
approach provides the temporal and spatial rapidity evolution of the typical hard momentum 
of the plasma partons, phard, and the plasma anisotropy, £. In a previous paper one of us [56] 
computed the thermal suppression of the T(ls) and Xbi states at LHC energies by folding 
together the aHydro temporal evolution of Ref. with results obtained in Ref. [IB] for 
the real and imaginary parts of the binding energy. In this paper, we extend this study to 
compute the suppression of T(ls), T(2s), T(3s), Xbi, and Xb2 states at both RHIC and LHC 
energies. 

The structure of the paper is as follows. In Section [2] we introduce the model potential 
we will use in order to compute the real and imaginary parts of the binding energies of the 
states under consideration. The potential used herein is an improved version of the one used 
in Refs. [56J and [16] and includes the effects of running coupling and an improved parame- 
terization of the numerical results for the short-range anisotropic potential. In Section [3] we 
briefly review the numerical method used to solve the Schrodinger equation. In Section [4] we 
review the aHydro dynamical model we use and discuss the qualitative features we expect 
to emerge based on the resulting dynamical evolution. In Section [5] we present the initial 
conditions we will use which include Glauber (or participant) scaling and a two-component 
model in which we use a linear combination of participant and binary collision scaling. In 
Section [6] we describe how we compute the nuclear modification factor Raa from the spatial 
and proper-time dependence of the real and imaginary parts of the binding energy. In Sec- 
tion [7] we detail how one can include the effect of feed down from excited states to compute 
the inclusive or "full" nuclear modification factor for the T(ls). In Section M we present 
our final results as a function of centrality, rapidity, and transverse momentum. Finally, in 
Section [9] we present our conclusions and outlook for future work. 

2. Setup and Model Potential 

In this section we specify the two potential models we consider in this paper. We consider 
the general case of a quark-gluon plasma which is anisotropic in momentum space. In the 
limit that the plasma is isotropic, the real part of the potentials used here reduces to the 
potential originally introduced by Karsch, Mehr, and Satz (KMS) [9J with or without an 
additional entropy contribution [12] and the imaginary part reduces to the result originally 
obtained by Laine et al [SO]- To begin the discussion we first introduce our ansatz for the 
one-particle distribution function subject to a momentum-space anisotropy. 
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2.1. Momentum- space anisotropic plasma 

The phase-space distribution of gluons in the local rest frame is assumed to be given by 
the following ansatz [4H 



/(t,X,p) = /iso (VP 2 +^(P • n) 2 /phard) 



where /i SO is an isotropic distribution which in thermal equilibrium is given by a Bose-Einstein 
distribution, £ is the momentum-space anisotropy parameter, and phard is a momentum scale 
which specifies the typical momentum of the particles in the plasma and can be identified 
with the temperature in the limit of thermal isotropic (£ = 0) equilibrium]^] The two param- 
eters Phard and £ can, in general, depend on proper time and position; however, we do not 
indicate this explicitly for compactness of the notation. The ansatz above is the simplest 
ansatz which allows for the breaking of symmetry in the Pt~Pl plane while maintaining local 
azimuthal symmetry in the transverse directions in momentum space. Note that one can use 
the same distribution to describe the quarks in the system [5TH5T?] and the quark self-energy 
in this case has been computed explicitly [60J. For our purpose, we are primarily interested 
in the gluon distribution since this will enter into the determination of the heavy quark 
potential; however, in the section on dynamics we implicitly assume the same distribution 
for quark degrees of freedom. 

Such a breaking of symmetry in the Pt~Pl plane arises naturally in a heavy-ion collision 
due to the rapid longitudinal expansion of the matter along the beamline direction and the 
parameter £ quantifies the degree of momentum-space anisotropy, 

2 [Pi) 

where p z = p • n and p^ = p — n(p • n) denote the particle momenta along and perpendic- 
ular to the direction n of anisotropy, respectively. For heavy ion collisions the anisotropy 
vector, n, lies along the beamline direction which we generally choose to lie along the z-axis. 

The energy-momentum tensor T^it, x, p) = (27r)~ 3 J d 3 p/p° p^p" f(t, x, p) for the dis- 
tribution function ([I]) is diagonal in the comoving frame and its components are [511 EI] 



1 The only place that we will assume thermal equilibrium herein is in the value of the isotropic Debye mass 



used in the heavy quark potential in Section 2.2.6 In principle, one could use another isotropic distribution 



function, in which case one would need to recompute the isotropic Debye mass. 
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£(Phard,0 =T TT = i + ^iso(Phard) , (3a) 

= 72.(0 ^tso(Phard) , 

TMPhard, = £ (T* + I™) = | + P iso ( Phard ) , (3b) 

= ^(OPiBoCPhard) , 

V L {p hard ,0 = ~T< = - I I Piso(Phard) , (3c) 

= 72 L (0^i so (Phard) , 

where Pi SO (phard) is the isotropic pressure and £i SO (Phard) is the isotropic energy density. In 
everything that follows we will use a conformal equation of state for which £ iso = 3Pi SQ . 

2.2. Model potential 

In this subsection we first review the derivation of the short range screened heavy-quark 
potential in the presence of finite momentum-space anisotropy. The full complex potential 
for an isotropic plasma was first obtained in Refs. [30j EI] • The calculation of the real part 
of the potential at finite anisotropy was first obtained in Ref. [12] and was later extended 
to include the imaginary part in Refs. 03TH5] . After this brief review we construct an 
analytic approximation to the real part of the heavy quark potential which allows us to 
compute the potential efficiently without having to resort to complicated two-dimensional 
numerical integration. As we will show, the resulting analytic approximation for the real 
part can be cast into the form of a Debye-screened Coulomb potential with a Debye mass 
which depends on the relative angle of the line connecting the quark and antiquark to the 
anisotropy direction. 

2.2.1. Integral expression for the real part of the short range potential 

One can determine the real part of the heavy-quark potential in the non-relativistic limit 
from the Fourier transform of the static gluon propagator. In an anisotropic plasma with a 
distribution function given by Eq. at leading order in the strong coupling constant one 
finds 02] 

V(v,0 = -g 2 C F / 7 ^ e ^A 00 (u; = 0,p,O, (4) 



2tt 



' ' / (27T) 3 ( p 2 + m 2 +m 2 )(p2 + m 2 ) _ m 4 > 
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where g is the strong coupling constant and Cp = (N 2 — 1)/ (2N C ) is the quadratic Casimir 
of the fundamental representation of SU(N C ). The masses in (|5| are given by 



2 D 

m n = 



' L ' " 2 arctanA/^ - arctan y ^P z _ j ) (5) 



ml = m 2 D 



(V? + (1 + Oarctanv^)(p 2 + £p 2 ) + £p, (p.VC + arctan 



2v^(l + 0(P 2 + ^l) 



(7) 



(2d 2 
P 2 1 + ^f . /7 , P,p 2 (2p 2 + 3^i) . Vt Pz 
7: 1=^ arctan Vc H ^ — arctan — . 



2 _ 7rm^P±|p| m 



+ p 
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with ttid being the isotropic Debye mass 



™ 2 D =-£ r^^, (io) 



27T 2 Jo C?p 

and p 2 = p 2 = p 2 ± + p 2 . The above expressions apply when n = (0, 0, 1) points along the 
z-axis; in the general case, p z and p_i_ get replaced by p ■ n and p — n(p • n), respectively. 
One can factorize the denominator of ^ by introducing 



2m 2 = M 2 ± yj M 4 — 4(m 2 3 (m 2 f + mty — mf) , (11) 
with M 2 = m 2 + m 2 + m 2 (57] . This allows us to write 

f d 3 v P 2 + ml + ml 

V(r, = -9 2 C F / ^ e'P' r . ^ ° . (12) 

7 (27r)-* (p 2 + m+)(p 2 + ml) 



In general one must evaluate (12) numerically. The integration can be reduced to a two- 
dimensional integral over a polar angle, 9, and the length of the three-momentum, p. How- 
ever, there can be poles in the integration domain due to the fact that m 2 _ can be negative for 
certain polar angles and momenta [SZjj^] These poles are first order and can dealt with using 
a principle-part prescription, however, evaluating this integral with the necessary precision 
requires on the order of 0.5 to 1 seconds per point. This presents a fundamental problem if 
one intends to evaluate the potential when solving the Schrodinger equation on large spatial 
lattices with on the order of 512 3 points. We are, therefore, motivated to find an efficient 
parametrization of the resulting potential based on a finite set of numerical evaluations. In 
order to do so, it is necessary to first consider various asymptotic limits of the potential. 



2 This is related to the presence of unstable collective modes in momentum-space anisotropic plasmas. 
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2.2.2. Asymptotic limits of the real part of the short range potential 

When £ = then m@ = m + = mo and all other mass scales are zero. As a consequence, 
we recover the isotropic Debye-screened Coulomb potential 

hmV(r,0 = VUr) = ~9 2 C F [ f P , f*\ = e~ f , (13) 



' 1> " J (2tt) 3 p 2 + m D Arrr 

where f = rmo- 

In the limit r — > for arbitrary £ one finds that the potential reduces to the vacuum 
Coulomb potential [42] 

(^V=-i^- (14) 

The same potential emerges for extreme anisotropy since all rrii — >■ as £ — >■ oo: 

lim^(r,0 = Kac(r) . (15) 

This is due to the fact that at £ = oo the phase space density /(p) from Eq. Q has 
support only in a two-dimensional plane orthogonal to the direction n of anisotropy. As a 
consequence, the density of the medium vanishes in this limit. 

2.2.3. Subleading terms in the small £ limit 

Having discussed the leading terms in the limits show above, we now discuss the sub- 
leading terms in the small £ limit. In the limit of small £ one finds that [57J 

mi = l + -(3cos20-l) , 
6 



m 2 = m 2 + m 2 = — cos 2$ , (16) 
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where rh = rn/rno and 9 is the angle with respect to the anisotropy vector n. As a result, 
one finds that 

lim V(t, = -g 2 C F I ^ J^- 2 , (17) 

where v = mc[l + | (3 cos 29 — 1)]. Expanding the integrand to leading order in £ and 
evaluating the resulting integrals one finds [42] 

lim V(r, = Mso(r) [1 - £F(f , 0)] , (18) 



where Vi SO (r) is the Debye-screened Coulomb potential in an isotropic medium (13), and the 
function J-(r, 9) = fo(f) + /i(f) cos(2#) with 



6(l-e r "')+f[6-f(f-3)] r f 2 
/oW = ^ = -6~48 + ---' (19) 



; = 6(l-eQ+r[6 + r(r + 3)] _ f 2 | 

4f 2 16 
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We can now define a ^-dependent screening mass in an anisotropic medium as the inverse 
of the distance scale r me( ±(9) over which |rV"(r)| drops by a factor of e: 

1 u^L (2i) 

To leading order in £ this leads to f me d = 1 — £^ r (r me d, #)• An approximate solution to this 
last equation gives [42J 

a 3 + cos 20 . , 

lim — ~ 1 - £ , 22 

where /i = r~l d . 

With this in hand we have an analytic approximation for the potential in the limit 
of small £, namely, that it is approximately a Debye-screened Coulomb potential with a 



^-dependent screening mass given by Eq. (22) such that 



hm^(r,£) ~ V iso (r) = ^ , ( 23 ) 

2.2.4- Subleading terms in the large £ limit 

We now turn our attention to the limit of large £. For general £ one can show that in an 
anisotropic plasma with a distribution function given by Eq. Q the particle number density 
can be factorized using a simple change of variables 

™(Phard, £) = J J^ffo X ' P) = ^/YZ^ ' ( 24 ) 

where n iso is the number density that would be obtained using the isotropic distribution 
function used in Eq. (jlj. Since in an isotropic system one can estimate the screening mass 
via m 2 D ~ n/T, we expect that in the large-£ limit one can will obtain /i 2 oc n(phard, £)/phard 
for the anisotropic screening mass, which leads to /i ~ £ -1 ' 4 J71.d in the large-£ limit. To see 
how this emerges analytically we return to the defining equation for the potential given in 
Eq. ([5]). In the limit of large £ one findsj^] 

to V M = V mV ) - -L-^J> J ^ . (25) 

We can compare this to the small screening-mass expansion of the isotropic potential Debye- 
screened Coulomb potential 

/d 3 p e ipr 

From the comparison we see that the anisotropic case can be obtained if we identify 



lim JL ~ xV 1/4 • (27) 
5^00 mo 2 



3 Note that the second integral below is infrared divergent and needs to be regulated; however, since we 
will only compare the coefficients of such integrals, we do not need to specify how it is regulated as long as 
we regulate it in the same manner in each case. 
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Figure 1: Comparison of the real part of the short range potential obtained from the analytic model 
specified in Eq. (28) and via direct numerical integration of Eq. (12). Panels (a)-(f) show the potential 
for different values of the anisotropy parameter as indicated in the lower left corner of each panel. In each 
panel the potential has been scaled by the vacuum Coulomb potential. Note that the vertical scale changes 
between panels. 



2.2.5. Model for the real part of the short range potential 

With both the small- and large-£ limits of the anisotropic screening mass in hand we 
can construct a model for the real part and then compare to direct numerical evaluation of 



the potential via Eq. (12). We find that the following form works well to reproduce the r 



dependence of the potential for all £. 



2 



i+e a 



+ + 



1/8 



with a = 16/tH, b = 1/2, d = 3/2, e 



(3 + £) 6 
1/3, and 



1 + 



c(0)(l + Q c 
(1 + ef) 



(28) 



c{6) 



3vr 2 cos(2#) + (9 + 4^3 - 4y/6) vr 2 + 64 (V6 - 3) 
4 (y/3 (y/2 - l) vr 2 - 16 (V6 - 3)) 



(29) 



The value of the parameter a in (28) guarantees that the large-£ form for the anisotropic 



screening mass (27) is reproduced. The expression for c(6) is determined by requiring that 
the small-£ limit (22) is reproduced. The coefficients b, d, and e were fit by hand in order 



to optimally reproduce the anisotropic short-range potential obtained by direct numerical 
integra 
ji/m D 



integration. In addition to reproducing these limits, the form (28) also guarantees that 
" 1//4 in the infinitely prolate limit of £ — > 



;i+e)- 



-1. We emphasize that the form 
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Figure 2: Comparison of the real part of the short range potential obtained from the analytic model 



specified in Eq. (28) and via direct numerical integration of Eq. (12) for £ = 1 



(28) is only a parametrization of the numerical results which is constructed in such a way 
as to guarantee the necessary asymptotic limits and to efficiently reproduce the potential 
obtained via direct numerical evaluation in an efficient manner. With this parametrization 
of \i in hand we can construct a model for the real part of the short range potential for all 

W{r)] = -^e-r, (30) 



with n given by Eq. ( 28 ) . 

In Fig. [I] we compare the model specified by Eq. (|30j) with results obtained by direct 
numerical integration for £ G {0.1, 1, 2, 10, 100, 1000} by plotting the ratio of the potential 
over the vacuum Coulomb potential. This is a very sensitive test of whether or not the 
parametrization is a good one and as can been seen from Fig. [I] works well over a very large 
range of possible plasma anisotropies. To see what the actual unsealed potential looks like 
in Fig. [2] we show the potential again; however, this time, we do not scale by the vacuum 



Coulomb potential. As can be seen from this figure, the model specified by Eq. (30) works 
extremely well allowing us to express the short-range anisotropic quarkonium potential as a 
Debye-screened Coulomb potential with an anisotropic screening mass [i. In the following 
subsection we will discuss the fact that one needs to model the long-range potential and 
construct a model for the potential at all scales. 

2.2.6. Model for the real part of the potential at all scales 

In order to make a realistic phenomenological model for quarkonium states one must 
consistently describe both short and long distance scales. Since heavy quark states are 
dominated by short distance physics at zero temperature they can be treated using heavy 
quark effective theory; however, as the temperature increases one expects the size of the 
states to increase causing the states to become sensitive to the long range part of the 
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potential. At zero temperature, since the velocity of the quarks in the bound state is 
small, quarkonium can be understood in terms of non-relativistic potential models such as 
the Cornell potential which can be derived directly from QCD using effective field theory 
P^0HT2] . A finite temperature extension of the Cornell potential might be provided by the 
KMS model [9] which describes the free energy of a static heavy quark-antiquark pair in an 
isotropic plasma via 



2/~( 

F(r,T) = -l±Z e -™ D r + e -m D rl /3-Q 

Airr m D L J 

where g is the strong coupling constant, o is the string tension, and mo is the isotropic 



Debye screening mass. Eq. (31) is a model for the action of a Wilson loop of size 1/T and 
r in the temporal and spatial directions, respectively (see [U2] and references therein). In 
the interest of spanning the possibilities for the real part of the potential we define potential 



model A by equating the real part of the potential with the free energy given in Eq. (31). 



However, in the general anisotropic case we must replace the isotropic screening mass by 



the anisotropic screening mass (28) to obtain 

SR[V A ] = F = -V" r + - [1 - e-» r ] , (32) 

where we have replaced g 2 CV/47r by a phenomenological parameter a in the screened coulomb 
contribution which will be adjusted to match lattice data. Here we take a = 0.385 which is 
consistent with the short range part of the heavy quark potential measured on the lattice 
[63] . For the isotropic Debye mass, m D , we use m 2 D = (1.4) 2 • N c (l + Nj/6) 47ra, j p 2 iard /3. 
The isotropic leading-order Debye mass is adjusted by a factor of (1.4) 2 in order to take into 
account higher-order corrections which have been measured in lattice simulations [64J. In 
the isotropic Debye mass we use a three-loop running for a s (65] with Aj^ = 344 MeV which 
gives a s (5 GeV) = 0.2034 in accordance with recent high precision lattice measurements of 
the running coupling constant [66J. For the scale of the running coupling we use 2ttT which 
is consistent with hard thermal loop calculations of quark-gluon plasma thermodynamics 
[6TI [68]. Finally, for the string tension we use a value of a = 0.223 GeV 2 which is again 
obtained from fits to lattice data [63]. In all cases we use N c = 3 since we are modeling QCD 
and take the number of contributing light quark flavors to be Nf = 2, which is appropriate 
for the temperature range considered herein^] 

As potential model B we will use the internal energy, U, of the states which has an 
entropy contribution added to it. To achieve this we calculate the full entropy S = —dF/dT 



using (32) and add T times this to the free energy (32), which leads to the internal energy 



U = F + TS. This procedure gives model B for the real part of the heavy quark potential 
R[V B ] = U = F-T^, (33) 

= --(l + /ir)e-' tr + — \l-e-» r ] -are-*"", (34) 
r /i 



4 If one uses instead Nf = 3 the isotropic Debye mass increases by ~ 6% which has only a small effect on 
the final results. 
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with jj, given by Eq. (28). In potential model B, we use the same parameters and Debye 
mass prescription as used in potential model A. 



2.2.7. Model for the imaginary part of the potential 

The imaginary part of the potential ^s[V] is obtained from a leading order perturbative 
calculation which was performed in the small anisotropy limit [S] . The resulting imaginary 
part of the potential is 



Q[V] = -a s C F T 



(r)-aMr,0) + Mr,0)) 



(35) 



where f = m D r, a s = # 2 /(47r), C F = (iV c 2 - 1)/(2N C ), and 



0(f) 
Vi(r,0) 



dz- 



z 



o 

oc 



dz- 



z 2 + iy 

z 



smizr 



z r 



1 



dz- 



o 



z 2 + 1 



2 

1 - 3 



sin 2 fl — ^ ; 1 n Q — 2 ' 



zr 



+ (1 -3 cos 2 6)G(r,z) 



(36) 
(37) 



cos 2 ^ 



zr 



+ (1 - 3 cos 2 9)G(f, z) 



with 9 being the angle from the beam direction and 

_ /A , fz cos(rz) — sin(f^) 
G(r,2) = 



\rz) 



(3f 



(39) 



For numerical efficiency three separate analytic expressions for S[V] which are valid in the 
small, medium, and large distance limits were determined and used in a piecewise fashion 
in their respective radii of convergence. 

2.2.8. Final Potential Models 

As mentioned above, here we consider two potential models, A and B, in which we identify 
the potential as coming from the free energy or internal energy, respectively. From both 
models discussed above we will additionally subtract a temperature- and spin-independent 
finite quark mass correction taken from Ref. [69J which improves the description of charm 
quark states at low temperatures, but is a small correction for bottom quarks. The final 
result for potential model A is 



V A = U[V A ]+i^[V] 



0.8 a 

2 ' 



Model A 



(40) 



with 9ft [Va] given by Eq. (32) and 5s\V] given by Eq. (35). The final result for potential 
model B is 

V B = sR[Vb] + &\y] - ^ , Model B (41) 



2 ' 

m 2 Q r 
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with 9ft [Vb] given by Eq. (34) and Q[V] given by Eq. (35). We note that both 9ft [Va] and 
9ft[Ve] reduce to the Cornell potential at T = and the short range part (r <C l/m D and 
r <C °f both reduces to the Coulomb potential, V = —a/r, at all temperatures, with 

a constrained by lattice data [63] . 



3. Solving the 3d Schrodinger Equation 

To solve the resulting Schrodinger equation we use the finite difference time domain 
method [701 E3] extended to the case of a complex- valued potential [IB]- Here we briefly 
review the technique. To determine the wave functions of bound quarkonium states, we 
must solve the time-independent Schrodinger equation 



E, 



H 



2m 



R 



+ V(x) + mi+m 2 



(42) 



on a three-dimensional lattice in coordinate space with the potential given by V = 9ft [V] + 



z£$[V] where the real and imaginary parts are specified in either Eqs. (40) and (41), respec- 
tively. Here, mi and m 2 are the masses of the two heavy quarks and m# is the reduced 
mass, mn = m\m2j {m\ + m-i). The index v on the eigenfunctions, <p v , and energies, E v , 
represents a list of all relevant quantum numbers, such as n, I, and m for a radial Coulomb 
potential. Due to the anisotropic screening mass, the wave functions are no longer radially 
symmetric if £ ^ 0. Nevertheless we still label the states as 15 (ground state) and IP (first 
p-wave excited state), respectively. 

To obtain the time-independent eigenfunctions we start with the time-dependent Schrodinger 
equation 

i^(x,t) = Hi/;(x,t), (43) 
which can be solved by expanding in terms of the eigenfunctions, <p v : 



V>(x,t) = X)' 



x e 



-iE v t 



(44) 



If one is only interested in the lowest energy states (ground state and first few excited states) 



an efficient way to proceed is to transform (43) and (44) to Euclidean time using a Wick 
rotation, r = it: 

%(x,r) = -^(x,r), (45) 



Or 



and 



^(x,r) = ^QA,(x)e-^ T • (46) 

V 

For details of the discretizations used etc. we refer the reader to Refs. (70J ITT] . 
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3.1. Finding the ground state 

By definition, the ground state is the state with the lowest energy eigenvalue, Eq. There- 



fore, at late imaginary time the sum over eigenfunctions (46) is dominated by the ground 
state eigenfunction 

lim ^(x, r) -> c o 0o(x)e- BoT . (47) 
Due to this, one can obtain the ground state wavefunction, <f) , and energy, E , by solving 



Eq. (45) starting from a random three-dimensional wavefunction, VWtiai( x > 0), an d evolving 
forward in imaginary time. The initial wavefunction should have a nonzero overlap with 
all eigenfunctions of the Hamiltonian; however, due to the damping of higher-energy eigen- 
functions at sufficiently late imaginary times we are left with only the ground state, 0o( x )- 
Once the ground state wavefunction (or any other wavefunction) is found, we can compute 
its energy eigenvalue via 

p / , \ {4>v\H\M fd s xtf,H<P v . . 

E - (T ■* 0O) = imr = !**<fyk ■ (48) 

To obtain the binding energy of a state, -E^bmd, we subtract the quark masses and the 
real part of the potential at infinity 

F ~ (F m m (M^EM] (AO) 



where 



V^O) = lim M[V(9,r)}, (50) 

|r|— >oo 



which is a purely real quantity. For an isotropic potential Voo is independent of the quantum 
numbers v and equal to either a /mo or 2a /mo for potential models A and B, respectively. 
In the anisotropic case, however, this is no longer true since the operator Voo(^) carries 
angular dependence. Its expectation value is, of course, independent of 6 but does depend 
on the anisotropy parameter £. 

3.2. Finding the excited states 

The basic method for finding excited states is to first evolve the initially random wave- 
function to large imaginary times, find the ground state wavefunction, (j) , and then project 
this state out from the initial wavefunction and re-evolve the partial-differential equation in 
imaginary time. However, there are (at least) two more efficient ways to accomplish this. 
The first is to record snapshots of the 3d wavefunction at a specified interval r snaps hot during 
a single evolution in r. After having obtained the ground state wavefunction, one can go 
back and extract the excited states by projecting out the ground state wavefunction from 
the recorded snapshots of ^(x, r) [701 fTTj . 

An alternative way to select different excited states is to impose a symmetry condition 
on the initially random wavefunction which cannot be broken by the Hamiltonian evolution 
|71j . For example, one can select the first p-wave excited state of the (anisotropic) potential 
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Figure 3: Real and imaginary parts of the T(ls), T(2s), and T(3s) binding energies as a function of the 
hard momentum scale, phard- The left panels show results obtained with potential model A (40 1 and the 



right panels show results from potential model B (41) 
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Figure 4: Real and imaginary parts of the Xbi and \b2 binding energies as a function of the hard momentum 



scale, phard- The left panels shows results obtained with potential model A (40) and the right panels show 



results from potential model B (41 1 



by anti-symmetrizing the initial wavefunction around either the x, y, or z axes. In the 
anisotropic case this trick can be used to separate the different excited state polarizations in 
the quarkonium system and to determine their energy eigenvalues with high precision. This 
high precision allows one to more accurately determine the splitting between polarization 
states which are otherwise degenerate in the isotropic Debye-screened Coulomb potential. 
Whichever method is used, once the wavefunction of an excited state has been determined, 



one can again use the general formulas (48) and (49) to determine the state's binding energy. 



3. 3. Results for the Binding Energies of Bottomonium States 

In Figs. [3] and [4] we show the real and imaginary parts of the T(ls), T(2s), T(3s), Xbh 
and Xb2 binding energies as a function of the hard momentum scale, phard, for £ G {0, 1, 20}. 



The left panels show results obtained with potential model A (|40|) and the right panels show 

t values of 
and T(3s) states 



results from potential model B (41). In each case we show three different values of £. For 

T(2s) 



the bottom quark mass we used m& = 4.7 GeV. For the T(ls) 
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£=0 


£=i 


State 


Potential A 


Potential B 


Potential A 


Potential B 


T(la) 


298 MeV 


593 MeV 


373 MeV 


735 MeV 


T(2s) 


< 192 MeV 


228 MeV 


< 192 MeV 


290 MeV 


T(3s) 


< 192 MeV 


< 192 MeV 


< 192 MeV 


< 192 MeV 


Xbi 


< 192 MeV 


265 MeV 


< 192 MeV 


351 MeV 


Xb2 


< 192 MeV 


< 192 MeV 


< 192 MeV 


213 MeV 



Table 1: Isotropic and anisotropic dissociation scales for the T(ls), T(2s), T(3s), Xbi, an d Xbi- Dissociation 
values were determined by finding the value of Phard when the real and imaginary parts of the state's binding 
energy become equal. 

we used a lattice size of N 3 = 256 3 with a lattice spacing of a = 0.125 GeV -1 ~ 0.025 fm 
giving a lattice size of L = Na ~ 6.3 fm. For the Xbi and Xb2 states we used a lattice size 
of TV 3 = 256 3 with a lattice spacing of a = 0.15 GeV -1 ~ 0.03 fm giving a lattice size of 
L = Na ~ 7.6 fm. Note that the fluctuations seen in some of the data points occur at 
values of Phard where the state is unbound. These fluctuations are due to poor convergence 
of the Schrodinger equation algorithm for unbound states. However, such fluctuations do 
not enter into our final results because, when the states are unbound (have a negative real 
part of their binding energy), then we use a large fixed decay rate for these states. Details 
of the precise prescription will be provided in Section |6j 

Defining the disassociation scale as the value of Phard at which the real and imaginary 
parts of the binding energy become equal, one finds the values listed in Table [TJ As can be 
seen from the figures and table one finds that the dissociation scale increases with increasing 
£ such that bottomonium states persist longer in a momentum-space anisotropic plasma. 
Binding energy data such as those presented in Figs. [3] and [4] will be used as input to our 
suppression calculation. 

4. Dynamical Model 

In order to describe the space-time evolution of the system we use "anisotropic hydrody- 
namics" (aHydro) which extends traditional viscous hydrodynamical treatments to cases 
in which the local momentum-space anisotropy of the plasma can be large EH]- The re- 
sult is a dynamical framework that reduces to 2nd order viscous hydrodynamics for weakly 
anisotropic plasmas, but can better describe highly anisotropic plasmas. In this paper we 
ignore the transverse expansion of the matter and model the system as a collection of de- 
coupled (1+1) -dimensional systems with different initial temperatures; however, we allow 
for the breaking of boost invariance. For such effectively one- dimensional dynamics which 
is homogeneous in the transverse directions, the aHydro approach provides the temporal 
and spatial rapidity evolution of the typical hard momentum of the plasma partons, Phard, 
the plasma anisotropy, £, and the four-velocity of the rest frame via a hyperbolic angle 

We briefly state the setup and final results of Ref. [55] for completeness. The starting 
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Figure 5: Dynamical parameters as a function of spatial rapidity using a strong coupling value of Ani]/S = 1. 
Shown are phard (left) and £ (right) with initial conditions £(to, ?) =0 and Phard( T o, ? = 0) = 540 MeV with 
t = 0.3 fm/c . The initial Phard rapidity dependence is given by a Gaussian profile specified in Eq. (62). 
Profiles at proper times r € {0.3, 2.1, 3.9, 5.7} fm/c are shown. 



point for the dynamical equations is to assume the same ansatz (fl| for the momentum-space 
anisotropic distribution distribution as was used to compute the heavy quark potential in 
the previous section. In the local rest frame of the plasma the ansatz has two parameters 
Phard and £. In the boost invariant case Phard and £ would be functions only of proper 
time; however, in the case of broken boost invariance both Phard and £ becomes functions of 
proper time, r, and spatial rapidity, q. The necessary dynamical equations can be obtained 
by taking moments of the Boltzmann equation in the relaxation time approximation [55] . 
The breaking of boost invariance requires that, in addition to phard and £, one must also 
specify the hyperbolic angle of the local rest frame of the flow. This can be accomplished by 
introducing two four-vectors, one of which specifies the four velocity of the local rest frame 
in lab frame, w M , and an additional four-vector, v^, which is orthogonal to u^, i.e. u^v^ = 0. 
This can be accomplished by introducing a hyperbolic angle d such that 

u* 1 = (coshtf(T,<;),0,0,sinhtf(T,^)) , (51a) 
u" = (sinhtf(r,<r), 0,0, cosher,?)) , (51b) 

where $(r, q) is the hyperbolic angle associated with the velocity of the local rest frame as 
measured in the lab frame [S3] . If the system were exactly boost-invariant then we would 
have i?(t, q) — q at all times. 

4-1. Moments of the Boltzmann Equation 

In order to obtain the necessary dynamical equations for Phard and £ we follow [HS] 
and take moments of the Boltzmann equation. For non-boost-invariant (l+l)-dimensional 
dynamics it suffices to take the zeroth and first moments and project the first moment with 
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Figure 6: Dynamical parameters as a function of spatial rapidity using a strong coupling value of 4ni]/S = 1. 
Shown are Phard (left) and £ (right) with initial conditions £(t , <t) =0 and Phard( r 0j ? = 0) = 350 MeV with 
To = 0.3 fm/c . The initial Phard rapidity dependence is given by a Gaussian profile specified in Eq. (62). 
Profiles at proper times r € {0.3,2.1,3.9,5.7} fm/c are shown. 



either w M or v^. The result is three coupled partial differential equations which give the 
proper-time and spatial-rapidity evolution of phard, £> and i?: 



^ _ 2(1 + ^\ _ _^ 9rPhard = 2A 



r 



Phard 



i-7e 3 / 4 (0^i + e' 



9. 



Q £ _|_ 4 ^hard + tanh(fl - 9) / g(Q ^ + ^ j^ghgA 

tanh(tf - 9 T « + 4 + i (ff 04 + 4 ^ 



Phard 

-fa 



7- V^l(0 
+ 1) (U 



g(g , , tanh(i?- ? ) 

^(0 



(52a) 



(52b) 



(52c) 



where 7£(£) and TZl(0 are defined in Eqs. (3a) and (3c), respectively. Note that in the 



derivation of the above equations it was assumed that the system consists of a plasma of 
massless particles which results in a conformal equation of state, i.e. £ iso = ?>P mo . 



The relaxation rate A appearing in the first equation (52a) is fixed by requiring that the 



equations reduce to the evolution equations of second order viscous hydrodynamics in the 
limit of small £. Doing so gives 



A 



2T(r) 2ft 1 / 4 (£)p hard 



5r/ 



5r/ 



(53) 
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Figure 7: Dynamical parameters as a function of spatial rapidity using a strong coupling value of Airrj/S = 
10. Shown are phard (left) and £ (right) with initial conditions ^(t , ^) = and Phard^o, 1 ? = 0) = 540 MeV 
with t = 0.3 fm/c . The initial Phard rapidity dependence is given by a Gaussian profile specified in Eq. (62). 
Profiles at proper times r g {0.3, 2.1, 3.9, 5.7} fm/c are shown. 



where f/ = rj/S is the ratio of the plasma shear viscosity to entropy density and we have 
mapped the equilibrium temperature to Phard and £ by requiring that the anisotropic and 
isotropic energy densities are the same, i.e. £ a mso(Phard, = ^iso(^), which upon using 
Eq. Q gives T = ^(Ofw 

We note, importantly, that since the relaxation rate A is proportional to phard, one expects 
that the relaxation to isotropic equilibrium is slower in regions where phard is lower. In addi- 
tion, we see that the relaxation rate is inversely proportional to fj which tells us that when 
the shear viscosity is small we expect to see larger plasma momentum-space anisotropics 
developing. In order to illustrate the dependence on initial temperature, in Figs. [5] and 
[6] we show the evolution of phard and £ in the case of a strong coupling shear viscosity of 
fj = 1/47T for two different assumed initial central temperatures of 540 MeV and 350 MeV, 
respectively. As can be seen from these two figures, as the initial temperature decreases, one 
sees larger momentum-space anisotropy as expected from Eq. (53). In order to illustrate 
the dependence on the assumed value of fj in Fig. [7] we show the case of f\ = 10/47T with 
an initial central temperature of 540 MeV. Comparing Figs. [5] and Fig. [7] we see that there 
is a dramatic increase in the developed momentum-space anisotropy when changing fj from 
1/47T to 10/47T. The result of these two dependences will be that we will see less suppression 
of the bottomonium states when Phard is l° w or V is large. 



5. Initial Conditions 

In this section we specify the type of initial conditions we use. We study both RHIC 
and LHC energies, therefore in this section we will present the general formulae which can 
be used in both cases. In the results section we will specify the specific initial temperatures, 
collision energies, starting proper times, etc. that we use in each specific case. 
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5.1. Transverse Coordinate Dependence 

In this paper we will consider collisions of symmetric nuclei, each containing A nucleons. 
We will study both participant and binary collision type initial conditions [72] using a 
Woods-Saxon distribution for each nuclei's transverse profile |73j. For an individual nucleon 
we take the nucleon density to be 



n A (r) 



n 



1 + e (r-R)/d 



(54) 



where no = 0.17 fm -3 is the central nucleon density, R = (1.12A 1 / 3 — fm is 

the nuclear radius, and d = 0.54 fm is the "skin depth". The density is normalized such 
that liniA^oo / d 3 rn A (r) = A, where A is the total number of nucleons in the nucleus. The 
normalization condition fixes no to the value specified above. From the nucleon density we 
first construct the thickness function in the standard way by integrating over the longitudinal 
direction, i.e. 



Ta{x,v) 



dzn A (\/x 2 + y 2 + z 2 ) . 



(55) 



With this in hand we can construct the overlap density between two nuclei whose centers are 
separated by an impact parameter vector b which we choose to point along the x direction, 
i.e. b = bx. We choose to locate the origin of our coordinate system to lie halfway between 
the center of the two nuclei such that the overlap density can be written as 



n AB (x,y,b) = T A (x + b/2,y)T B (x - b/2,y) . 



(56) 



The overlap density will be used later as the probability weight for bottomonium production 
and our "two-component" initial condition. Another quantity of interest is the participant 
density which is given by 



n paxt (x,y, b) = T A (x + b/2,y) 



(J NN T B (x - b/2,y) 



B 



+ T B (x-b/2,y) 



a NN T A (x + 6/2, y) 
A 



■ (57) 



For LHC collisions at a/saT/v = 2.76 TeV we use a^N — 62 mb and for RHIC collisions at 
y/SNN = 200 GeV we use 0"^^ = 42 mb. From the participant density we construct our 
first possible initial condition for the transverse phard profile at central rapidity by taking 
the third root of the rescaled n part 



part 
Phard,0 



■n pa ,rt(x,y,b) 



nl/3 



(5f 



n part (0,0,0) 

where T is the central temperature obtained in a central collision between the two nuclei. 
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As an alternative initial condition for Phard 
one could use the number of binary collisions 
which is defined as 

n co \i(x, y, b) = a NN n AB (x, y, b) . (59) 

Comparisons with RHIC data show that it is 
necessary to add an admixture of n co a to the 
participant, or wounded- nucleon, scaling. We 
will consider such an admixture as our second 
possibility by defining 



n mix {x,y,b) 



: (1 - a)n paxt (x,y,b) 
+an co \\(x,y,b) , (60) 



with a = 0.145 as fit by the PHOBOS Collab- 
oration [73]. This gives a second possibility 
for the initial condition for Phard a t central ra- 
pidity 
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Figure 8: Comparison of initial transverse phard 
profile given by n part scaling and n co ii scaling. A 
value of (T^v n = 62 mb was used and we show the 
case of a central collision, i.e. 6 = 0. 



mix 
Pha,rd,0 



(x,y,b) 



n. 



1/3 



(61) 

(0,0,0). 

Note that T should be adjusted so that both initial conditions give the same particle density 
at central rapidity when integrated over the transverse plane. For a = 0.145 we find that 
T mix = 1.079 T part at LHC energies and T mix = 1.065 T part at RHIC energies. These values 
will be used in the results section when we discuss the initial condition dependence of our 
results. 



5.2. Spatial Rapidity Dependence 

In the previous subsection we fixed two possible prescriptions for the transverse tem- 
perature profile. Since we allow for the breaking of boost-invariance, we also need to give 
the spatial-rapidity dependence in order to complete our specification of the full three- 
dimensional initial temperature profile. For the number density profile in spatial rapidity 
(<j) we use a Gaussian which successfully describes experimentally observed pion rapidity 
spectra from AGS to RHIC energies [7S14T9] and extrapolate this result to LHC energies. 
The parametrization we use is 



n(q) = n exp[ 4 



2a? 



(62) 



with 



In {^s NN /2m p ) 



(63) 
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where c s is the sound velocity, m, p = 0.938 GeV is the proton mass, y/s^N is the nucleon- 
nucleon center-of-mass energy, and n is the number density at central rapidity. We have 
added a multiplicative factor of 0.64 to adjust for broadening of the distribution in rapidity 
as a function of proper time since the fits, e.g. from [79], were to the final state spectra 
rather than initial state spectra. In this paper we will use an ideal (conformal) equation of 
state for which c s = in natural units. 

5.3. Full Three- Dimensional Initial Conditions 

We can use Eq. (62) to determine the initial phard rapidity dependence by taking the 
third root of the number density. Putting this together with the two possibilities for the 



transverse temperature dependence determined in the Section 5.1 we can now specify the 
full three-dimensional initial temperature profile. Depending on whether we use the number 
of participants (n part ) or two component model (n m j x ) scaling we have two possible initial 
Phard profiles: 



Phard,0 ~~ ^0 



n 



part 



{x,y,b) e" ?2 /( 2<J ?) 



nl/3 



n 



part 



(0,0,0) 



; Initial Condition I , 



(64) 



II _ rp 

Phard,0 ~ J 



n mix {x,y,b)e « 2 /( 2(j2 ) 



71, 



,(0,0,0) 



1/3 



Initial Condition II 



(65) 



5.4- Allowing for initial momentum- space anisotropy 

If the initial momentum-space anisotropy is assumed to be zero, i.e. £o 



0, then Eqs. (64) 



and (65) can be used without modification. However, if £ 7^ one should require that the 
same initial density profile is obtained. Using the fact that n(phard ; = ^iso(Phard) / \/l + £ oc 
Phard/v / ^+T one nn ds that this requires T (£ ) = (1 + £,) 1/6 T - lso . 

We must note, however, for completeness sake, that one could also have a non-trivial 
dependence of the initial anisotropy on the transverse direction and spatial rapidity. In fact, 
one expects that towards the transverse and longitudinal edges of the plasma that the initial 
momentum-space anisotropies should be larger; however, at this point in time there is no 
first principles calculation of the xi and ? dependence of £ at the earliest times after the 
collision, so here we will choose the simplest possibility, which is that it is a constant and 
equal to zero. We will explore the possibility of finite initial momentum-space anisotropy in 
future works. 



6. Computing the suppression factor 

The aHydro time evolution gives us Phard and £ as a function of proper time, transverse 
coordinate xj_, and spatial rapidity Solution of the Schrodinger equation gives us the 
real and imaginary parts of the binding energy of a given state as a function of phard and 
£. Putting this together gives us the real and imaginary parts of the binding energy as a 
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function of proper time, transverse coordinate xj_, and spatial rapidity q: ^ft[E hind (r, xj_, <;)] 
and Q^-EbmdtVj xj_, <;)], respectively. 

If the real part of the binding energy is positive, then the state is bound. If the real part 
of the binding energy is negative, then the state is unbound. The imaginary part of the 
binding energy will give us information about the decay rate of the state in question. To 
see the exact relationship we can compute the quantum mechanical occupation number as 
a function of proper time 



n v (r) = (0*(r,x)0„(r,x)}, 

= ((^(x) e -^)*(^(x) e -^)), 

= (0:(x)0 1) (x))e 23 ^\ 

= n°e 23 ^, (66) 

where in the last line we have identified n° v = (0* (x)0 1) (x)). In order to connect this to the 
decay rate, T, we note that T is defined empirically through n v (t) = n° v exp(— Yr) so that 



we can identify Y = — 2$s[E\. Finally, from Eq. (49) we have S[i?bind] = — ^[E] so that 



r / A _\ 29[£ bind (r,x ± ,<)] sft[£ bind (r,x ± ,^)] > , . 

1 ' ±,q) " I 10 GeV R^foxj.,?)] < ^ 

The value of 10 GeV in the second case is chosen to be large in order to quickly suppress 
states which are fully unbound. We have checked the sensitivity of our results to this value 
and find that there is very little dependence on this number as long as it is greater than 1 
GeV such that the states are suppressed quickly within the plasma lifetime. In addition, we 
set the width to zero if the imaginary part of the binding energy is less than zero. Negative 
values of the imaginary part of the binding energy occur only at large values of £ and are 
a result of the small-£ expansion being applied outside of its range of applicability. Since 
large £ corresponds to a (nearly) free streaming plasma, one expects that the widths should 
return to their vacuum values (~ keV) justifying this choice. 

We can integrate the instantaneous decay rate, Y, over proper-time to extract the di- 
mensionless logarithmic suppression factor 

r T f 

C(Pr,x±,0 = 6(77 - r form (p T )) / drY(T,x ± ,q) , (68) 

Jmax(r form (p T ),T ) 

where rf or m(pr) is the lab-frame formation time of the state in question. The formation 
time of a state in its local rest frame can be estimated by the inverse of its vacuum binding 
energy [80] . In the lab frame the formation time depends on the transverse momentum of 
the state via the gamma factor Tf orm (p;r) = 7 r fo rm = ET T f OTm /M where M is the mass of the 
relevant state and T f ° rm is the formation time of the state in its local rest frame. For the 
formation times for the T(ls), T(2s), T(3s), Xbi an d Xb2 states we take r f ° rm = 0.2 fm/c, 
0.4 fm/c, 0.6 fm/c, 0.4 fm/c, and 0.6 fm/c, respectively. 

We take the initial proper time r for plasma evolution to be r = 0.3 fm/c at both RHIC 
and LHC energies. The final time, r/, is defined to be the proper time when the local energy 
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T(ls) Production 


Mechanism 


% ± Stat ± Sys [82J 


fi used herein 


Direct Production 


50.9 ± 8.2 ± 9.0 


0.51 


T(2s) decay 


10.7 ± 7.7 ± 4.8 


10.7 


T(3s) decay 


0.8 ± 0.6 ± 0.4 


0.8 


Xbi decay 


27.1 ± 6.9 ± 4.4 


27 


Xb2 decay 


10.5 ± 4.4 ± 1.4 


10.5 



Table 2: Feed down fractions extracted from experiment |82j including errors (middle column) and the value 
chosen for use herein (right column). Values of fi are constrained such that /« = 1- 



density becomes less than that of an N c = 3 and Nf = 2 ideal gas of quark and gluons 
with a temperature of T = 192 MeV. At this energy density, plasma screening effects are 
assumed to decrease rapidly due to the transition to the hadronic phase and the widths of 
the states will become approximately equal to their vacuum widths. 



From £ obtained via Eq. (68) we can directly compute the suppression factor Raa 

Raa(pt^±^) = e~^ T ^ . (69) 

For averaging over transverse momenta and implementing any cuts necessary we assume 
that all states have a 1/ E\ spectrum which is consistent with the high-pr spectra measured 
by CDF |81j. Integrating over transverse momentum given p-p-cuts pr,mm and pT.max we 
obtain the p^-cut suppression factor 



For implementing cuts in centrality we compute Raa for finite impact parameter b and map 
centrality to impact parameter in the standard manner. For the cuts over centrality and 
rapidity, we use a flat distribution. 

In order to compare with experimental observations we should finally average Raa( x -±,s) 
over Xj_. For this operation we use a production probability distribution which is set by the 



overlap density specified in Eq. (56) 



J x± dx ± n A A{x±)RAA{x±,<;) 

(Raa{<;)) = r~r - f — x • (71) 



7. Excited State Feed Down 

Since a certain fraction of T(ls) states produced in high energy collisions come from the 
decay of excited states, when computing the full (inclusive) Raa for the T(ls) one must also 
consider the suppression of the excited states which decay or "feed down" to it. In order to 
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fix the feed down fractions we use data from a/s =1.8 TeV pp collisions at CDF [H2] with 
a cut > 8.0 GeV/c. The resulting feed down fractions are listed in Table [2] 

Based on these numbers, we can construct the full (or inclusive) T(ls) Raa including 
the effect of the suppression of excited states via 

RfX[T{ls)\= /<^> ( 72 ) 

i £ states 

where R^aa is the direct suppression of the i th state and the production fractions, fi, are 
given in Table [2] 



8. Results for Raa 

In this section we present our main results which consist of the suppression factors Raa 
for the T(ls), T(2s), and T(3s), Xbi an d Xb2- We will present each state's suppression factor 
as a function of centrality (number of participants) and rapidity. We will then compute the 
inclusive Raa for the T(ls) including the feed effect as described in Section [7j To close the 
section we will present the inclusive Raa for the T(ls) as a function of transverse momentum 
and investigate the sensitivity to the choice of the type of initial conditions used. 

8.1. Suppression at RHIC Energies 

The highest energy RHIC runs collide gold nuclei at a collision energy of ^/snn = 200 
GeV. In this subsection we will focus on the resulting using wounded-nucleon (or participant) 
scaling for the initial condition with u^n — 42 mb. Fixing the initial time for the aHydro 
evolution to r = 0.3 fm/c and requiring that the final charged particle multiplicity is fixed to 
dNch/dy = 620, we find that for Anrj/S = {1, 2, 3} we must fix the initial central temperature 
for a central collision to be To = {442,433,428} MeV. The decrease of the initial central 
temperature with increasing rj/S is a result of the fact that one has more entropy generation 
as rj/S increases. As a result, it is necessary to lower the initial temperature in order to 
allow for particle production. 

In Fig. [9] we show the predicted suppression factor Raa for the T(ls), T(2s), T(3s), 
Xbu an d Xb2 states as a function of the number of participants (left) and rapidity (right). 



The top row uses potential model A (40) and the bottom row uses potential model B 



(41). In all plots we used ^/s^n = 200 GeV, assumed a shear viscosity to entropy density 



ratio of Atrrj/S = 1, and implemented cuts of < pr < 20 GeV and and (left) rapidity 



\y\ < 0.5 (right) centrality 0-100%. As can be seen from this figure, potential model A (40) 



provides much more suppression than potential model B (41), both as a function of number 
of participants and rapidity. In both cases we see clear signs of sequential suppression, with 
the higher excited states having stronger suppression than the ground state. However, we 
note that even for states that are melted at relatively low central temperatures, we still 
obtain a non-vanishing suppression factor for these states. This is due to the fact that near 
the edges, where the temperature is lower, one does not see suppression of the states. Upon 



performing the geometrical average prescribed in Eq. (71) we see that a large fraction of the 
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Figure 9: RHIC suppression factor Raa for the T(ls), T(2s), T(3s), Xbi, and Xbi states as a function of 
the number of participants (left) and rapidity (right). The top row uses potential model A (40 1 and the 
bottom row uses potential model B (41 1. In all plots we used ^/snn = 200 GeV, assumed a shear viscosity to 
entropy density ratio of Airq/S = 1, and implemented cuts of < pr < 20 GeV and (left) rapidity \y\ < 0.5 
(right) centrality 0-100%. 



states produced can survive even when the central temperature of the plasma is above their 
naive dissociation temperature. 



In Fig. 10 we show the inclusive suppression factor i?^[T(ls)] obtained using the feed 
down prescription presented in Section [7j As can be seen from these figures, potential 
model A (free energy) predicts much stronger suppression than potential model B (internal 
energy). As we can see the result has a significant dependence on the assumed shear viscosity 
to entropy density ratio. This could, in principle, be used to constrain r]/S from RHIC data 
on bottomonium suppression. 

8.1.1. Raa for T(ls + 2s + 3s) and comparison to STAR data 

Due to limited statistics and resolution the STAR Collaboration does not report separate 
suppression factors for the T(ls), T(2s), and T(3s) states. Instead, they compute an effec- 
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Figure 10: RHIC inclusive or "full" suppression factor Raa for the T(ls) including feed down effects. The 
three different lines correspond to different assumptions for the shear viscosity to entropy ratio Airrj/S £ 
{1, 2, 3}. In all plots we used y/sjyN = 200 GeV and implemented cuts of < pt < 20 GeV and and (left) 
rapidity \y\ < 0.5 (right) centrality 0-100%. 



tive total suppression of all three states by integrating the counts in a dielectron-invariant 
mass window which encompasses all three states 



m + An-*-, ^AA 



f m+ dm 



n: 



R AA [T(ls + 2s + 3s)] = , (73) 

where m_ and m + are the dielectron pair invariant masses which cover the T(ls), T(2s), 
and T(3s) spectral peaks, e.g. m_ = 8.5 GeV and m + = 11 GeV. If the spectral peaks have 
approximately the same width and are well separated, as is the case with these three states, 
then one finds that to good approximation 

Rjuffjl. + 2s + 3.)] * ^ |Y(1S)1 + g^ggfjl + ^[T(3 3) ] 

1 + c 2s + c 3s 
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Figure 11: RHIC T(ls + 2s + 3s) suppression factor determined via Eq. (74 1 compared with experimental 
data from the STAR Collaboration [S3] . The three different lines correspond to different assumptions for the 
shear viscosity to entropy ratio Airrj/S € {1, 2, 3}. In all plots we used ^snn = 200 GeV and implemented 
cuts of < p T < 20 GeV and \y\ < 0.5. 



where ci s and c^s are the ratios of the T(2s) and T(3s) states' background subtracted p-p 
peak heights to the T(ls) state's background subtracted p-p peak height, respectively. From 
preliminary LHCb results |84j in the dimuon channel one finds c<i s — 0.24 and c^ s ~ 0.11. 
These values are consistent with CMS measurements of the T(2s)/T(ls) and T(3s)/T(ls) 
cross section ratios [85] - We will use these values assuming that they are a good approxi- 
mation to the relative p-p peak heights in the dielectron channel. 



In Fig. 11 we plot Raa[Y(1s + 2s + 3s)] as determined using Eq. (74) and compare 
with experimental data from the STAR Collaboration (83] . As can be seen from this figure, 
potential model A (free energy) gives too much suppression when compared to RHIC data. 
One could argue that there could be some enhancement from regeneration; however, at 
RHIC, in particular, the number of bottom and anti-bottom quarks generated on an event- 
by-event basis is incredibly small and therefore regeneration due to recombination of the 
bottom and anti-bottom quarks is highly improbable. Potential model B, on the other 
hand, does a very good job of reproducing the existing STAR data for Raa[Y(1s + 2s + 3s)}. 
From the right panel we can obtain an estimate of r]/S: 0.08 < rj/S < 0.24. Unfortunately, 
a more accurate determination will require more data with reduced statistical errors. 

8.2. Suppression at LHC Energies 

The current LHC runs collide lead nuclei at a collision energy of y/s^M = 2.76 TeV. 
In this subsection we will focus on the resulting using wounded-nucleon (or participant) 
scaling for the initial condition with <jnn = 62 mb. Fixing the initial time for the aHydro 
evolution to tq = 0.3 fm/c and requiring that the final charged particle multiplicity is fixed 
to dN c h/dy = 1400, we find that for 4irr]/S = {1,2,3} we must fix the initial central 
temperature for a central collision to be T = {567, 550, 539} MeV. As before, the decrease 
of the initial central temperature with increasing rj/S is a result of the fact that one has 
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Figure 12: LHC suppression factor Raa for the T(ls), T(2s), T(3s), Xbii an d Xb2 states as a function of 
the number of participants (left) and rapidity (right). The top row uses potential model A ( |40| and the 
bottom row uses potential model B (41 ). In all plots we used ^Jsmn — 2.76 TeV, assumed a shear viscosity 
to entropy density ratio of Airrj/S = 1, and implemented cuts of < pr < 20 GcV and (left) rapidity 
\y\ < 2.4 (right) centrality 0-100%. 



more entropy generation as rj/S increases. 

In Fig. 12 we show the predicted suppression factor Raa for the T(ls), T(2s), T(3s), 
Xbi, and Xb2 states as a function of the number of participants (left) and rapidity (right). 
The top row uses potential model A (40) and the bottom row uses potential model B (41). 
In all plots we used ^snn = 2.76 TeV, assumed a shear viscosity to entropy density ratio 
of 4:711] /S = 1, and implemented cuts of < pr < 20 GeV and (left) rapidity \y\ < 2.4 
(right) centrality 0-100%. As can be seen from this figure, as was the case at RHIC energies, 
potential model A (40) provides much more suppression than potential model B (41), both 
as a function of number of participants and rapidity. In both cases we see clear signs 
of sequential suppression, with the higher excited states having stronger suppression than 
the ground state. However, we once again note that even for states that are melted at 
relatively low central temperatures, we still obtain a non-vanishing suppression factor for 
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Figure 13: LHC inclusive or "lull" suppression factor Raa for the T(ls) including feed down effects 
compared to experimental data are from the CMS Collaboration [53]. The three different lines correspond 
to different assumptions for the shear viscosity to entropy ratio Atttj/S E {1,2,3}. In all plots we used 
y/SNN = 2.76 TeV and implemented cuts of < pr < 20 GeV and (left) rapidity \y\ < 2.4 (right) centrality 
0-100%. 



these states. This is due to the fact that near the edges, where the temperature is lower, one 
does not see suppression of the states. Upon performing the geometrical average prescribed 
in Eq. ( 71 ) we see that a large fraction of the states produced can survive even when the 
central temperature of the plasma is above their naive dissociation temperature. 

In Fig. 13 we show the inclusive suppression factor i?^[T(ls)] obtained using the feed 
down prescription presented in Section [7j As can be seen from these figures, potential 
model A (free energy) predicts much stronger suppression than potential model B (internal 
energy). Comparing to the available CMS data [86] we see that, as was the case at RHIC 
energies, potential model B (internal energy) does a much better job of reproducing the 
data than potential model A (free energy) both as a function of centrality and rapidity. 
Using the potential model B results we can obtain an estimate for r)/S at LHC energies: 
0.08 < rj/S < 0.24 which is the same range obtained from the STAR data obtained with 
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Figure 14: LHC inclusive or "full" suppression factor Raa for the T(ls) including feed down effects 
as a function of transverse momentum compared to experimental data are from the CMS Collaboration 
|86| . The three different lines correspond to different assumptions for the shear viscosity to entropy ratio 
4nri/S <E {1,2,3}. For the plot we used ^/snn = 2.76 TeV and implemented cuts of < px < 20 GeV, 
\y\ < 2.4, and centrality 0-100%. 



gold-gold collisions at lower energies. As before, more precisely determining rj/S will require 
more data from the LHC which should be forthcoming in the near future. 

8.3. Transverse momentum dependence 

In Fig. [l4] we plot the minimum bias (centrality 0-100%) full suppression factor for 
the T(ls) including feed down effects as a function of transverse momentum. Since we 
ignore the transverse expansion of the matter created in the heavy ion collision, the only 
Pt dependence which is included comes from the formation time effect. One expects based 
on this that higher pt states will have weaker suppression since, in the lab frame, they are 
formed at a later proper-time when the plasma is cooler. This expectation is borne out by 



Fig. [T4J however, as can be seen from this figure there is only a weak p^-dependence of the 
result. This is to be contrasted with the relatively much larger px dependence of the CMS 
results. Looking forward we note that there are two additional places where a momentum 
dependence could enter the final results: (1) an intrinsic velocity dependence of the damping 
rate itself and (2) the effect of heavy quark states being nearly free streaming in the soft 
background. It has been shown [87] that adding a finite velocity relative to the medium 
affects the heavy quark potential, so this would indeed be something that one will need to 
investigate in future work and may help to improve agreement with the experimental data. 
The second effect will require the simultaneous solution of transport equations for nearly 
free streaming heavy quark states and the soft sector. 
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Figure 15: RHIC (left) and LHC (right) inclusive suppression factor Raa for the T(ls) including feed down 
effects compare to STAR [S3] and CMS [86J data. In both plots we have fixed Airrj/S = 2. Collision energies 
and cuts applied are indicated in each figure. The solid black line is the result obtained assuming wounded 
nucleon initial conditions and the dashed red line is the result obtained used a two component model with 
a = 0.145. 



8.4- Dependence on the choice of initial condition type 

As detailed in Section [5] we consider two types of initial conditions: (I) the wounded 



nucleon model specified in Eq. (64) and (II) a two- component model which consists of an 



admixture of participant and binary scaling specified in Eq. (65). In Fig. 15 we show the 
results obtained for Raa[^ (Ls + 2s + 3s)] at RHIC energies and the full (or inclusive) Raa for 
the T(ls). In both plots we have assumed Anrj/S = 2. Because changing the initial condition 
type affects particle multiplicities we have adjusted the initial temperature at RHIC energies 
from 433 MeV to 461 MeV and at LHC from 567 MeV to 612 MeV in order to keep the 
charged particle multiplicity fixed at dN ch /dy = 620 and dN ch /dy = 1400, respectively. As 



can be seen from Fig. 15 , for peripheral collisions there is a larger dependence on the choice 
of initial condition type, while for central collisions the result obtained is not much affected 
by the choice of initial condition. This is to be contrasted with the dependence of the result 
on the assumed value of r)/S which affects the suppression at all centralities. This leaves 
hope that one can disentangle the initial condition effect and the effect of the assumed value 

of Tj/S. 

9. Conclusions and Outlook 

In this paper we considered the suppression of bottomonium states in ultrarelativistic 
heavy ion collisions. We computed the suppression as a function of centrality, rapidity, and 
transverse momentum for the states T(ls), T(2s), T(3s), Xbi, and Xb2- Using this informa- 
tion, we then computed the inclusive T(ls) suppression as a function of centrality, rapidity, 
and transverse momentum including feed down effects. Calculations were performed for both 
RHIC = 200 GeV Au-Au collisions and LHC v /s^7 = 2.76 TeV Pb-Pb collisions. 
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Our calculations build upon a concerted theoretical effort to understand recently obtained 
RHIC and LHC data on bottomonium suppression [SU I8"8"H9~T] . 

We studied two different potential models which were based on the heavy quark free 
energy (A) and internal energy (B). We found that the potential based on the free energy 
gives too much suppression when compared to the available experimental data at both RHIC 
and LHC energies. On the other hand, results obtained from the potential model that was 
based on the internal energy seem to be in reasonably good agreement with data obtained 
at both collision energies. We are therefore led to conclude that one should not use potential 
models based on the free energy. From the comparison of our theoretical results obtained 
using the potential based on the internal energy and data available from the STAR and 
CMS Collaborations we were able to constrain the shear viscosity to entropy ratio to be in 
the range 0.08 < rj/S < 0.24. We find that our results are consistent with the creation of a 
high temperature quark-gluon plasma at both RHIC and LHC collision energies. 

That being said, it is worrisome that one sees such a strong dependence of the results 
on the potential model used. However, herein we find that at both RHIC and LHC energies 
a potential based on the internal energy seems to better describe the available data with 
values for the shear viscosity to entropy ratio which are consistent with those determined 
from bulk collective flow. The dependence on the potential used emphasizes the need for a 
concerted theoretical effort to better determine the heavy quark potential analytically via 
finite temperature effective field theory methods and/or numerically via lattice QCD studies. 
This will require determination of the both the real and imaginary parts of the potential 
at short and long distances and also the dependence on the momentum-space anisotropy of 
the plasma partons. The calculation of the short range part of the potential for arbitrary 
momentum-space anisotropy is currently underway. 

In future work we also plan to include the effect of allowing heavy quark states to have a 
flow which is decoupled from the soft medium and to include the effect of finite velocities on 
the heavy quark decay rate. This will include the addition of full 3+ld aHydro evolution 
so that we can simultaneously describe elliptic flow and bottomonium suppression. It would 
also be interesting to investigate the behavior of heavy quarkonium widths near T c using 
an AdS/QCD model. Finally, it will also be necessary to investigate the possibility of pair 
recombination due to residual spatial correlations among suppressed pairs [921 [93]. How 
these future investigations will affect the quoted range for rj/S is a critical open question 
which will need to be addressed. We leave these interesting questions for future work. 
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